Order 1/N corrections to the time-dependent Hartree approximation for a system of 

N + 1 oscillators 



Bogdan Mihailatl and John F. Dawsonl 

Department of Physics, University of New Hampshire, Durham, NH 03824 

Fred Coopeii 

Theoretical Division, MS B285, Los Alamos National Laboratory, Los Alamos, NM 87545 

(February 1, 2008) 

We solve numerically to order 1/N the time evolution of a quantum dynamical system 
of TV oscillators of mass m coupled quadratically to a massless dynamic variable. We use 
Schwinger's closed time path (CTP) formalism to derive the equations. We compare two 
methods which differ by terms of order 1/N 2 . The first method is a direct perturbation 
theory in 1/N using the path integral. The second solves exactly the theory defined by the 
effective action to order 1/N. We compare the results of both methods as a function of N. 
At N = 1, where we expect the expansion to be quite innacurate, we compare our results to 
an exact numerical solution of the Schrodinger equation. In this case we find that when the 
two methods disagree they also diverge from the exact answer. We also find at N = 1 that 
the 1/N corrected evolutions track the exact answer for the expectation values much longer 
than the mean field (N — oo) result. 
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I. INTRODUCTION 

The large N approximation has a long history in both 
statistical mechanics and quantum field theory, mostly in 
determining the phase structure of various theories jjj . It 
is only recently that this approximation has been used to 
study the dynamical evolution of various systems, and 
at present only the lowest order in the large N expan- 
sion has been considered. In leading order, the large 
N expansion is equivalent to using a gaussian density 
matrix, and therefore two particle scattering effects are 
included only indirectly. The leading order in the large 
N approximation is closely related to a time dependent 
Hartree approximation. The exact connection between 
these methods is discussed in detail in ref . j2| . Although 
interesting results in lowest order have been obtained for 
pair production from strong fields &5] as well as the 
evolution of a chiral phase transition [5], the important 
effects of the direct two particle elastic scattering, which 
determines the thermalization time scale of the plasma, 
is not included in the mean field approximation. In or- 
der to compare time scales for rethermalization with, say, 
the plasma oscillation frequency, as well as the expansion 
time of an evolving plasma produced in a heavy-ion col- 
lision, one needs to go to next order in the 1/N expan- 
sion. The property of the 1/N expansion relevant here 
is that connected 2n-point Green's functions first appear 
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with Gin of order 1/7V™ -1 . Thus direct scattering of two 
particles first occurs at order 1/N. At lowest order in 
1/N, the equations one has to solve are differential equa- 
tions. At next order one gets integro-differential equa- 
tions which depend on the time history of the system. 
This requires new numerical methods to ensure that the 
update conserves energy. 

There are two different ways one can determine the 
1/N corrections. The first method is to iterate the solu- 
tion of the lowest order calculation in a standard pertur- 
bative fashion. However, one might hope that it might 
be more accurate to first Legendre transform the action 
obtained to order 1/N and to obtain a new effective ac- 
tion, which differs by terms of order 1/N 2 from the first 
method. In the second method, one evolves directly the 
equations of motion obtained from the effective action. 
By having these two different methods, one has an upper 
bound on the accuracy of the 1 /N expansion. When the 
two methods diverge, this is a signal that 1/N 2 correc- 
tions are important. At iV = 1 this divergence is very 
close to the place where these two methods diverge from 
the exact answer. We find that the method based on the 
effective action is actually less stable, since solutions be- 
come unbounded earlier. However this occurs much after 
the method is unreliable. 

In quantum field theory it is not possible to compare 
the 1/N expansion with an exact calculation of a dynam- 
ical evolution because of the large number of degrees of 
freedom. Thus we thought it appropriate to study a sim- 
ple quantum mechanical example which we have studied 
before in the lowest order mean field approximation . 
The advantage of this simple model is that, at least for 
N = 1 , comparisons can be made with a direct numerical 
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simulation of the exact problem j?j . Unfortunately, even 
for the quantum mechanics case, going beyond N = 1 
is not numerically feasible for the exact problem. So we 
are testing our large N expansion in a regime where it is 
not expected to work very well. However, in spite of this 
short coming, we find that by adding the 1/N correction 
terms, our approximation tracks the exact answer for ex- 
pectation values a factor of two longer than the lowest 
order approximation — which is encouraging. 

As far as we know, this is the first attempt to use 
Schwinger's CTP formalism in a calculation which is not 
a perturbation expansion in the coupling constant. The 
methods of solving the resulting Volterra-like equations, 
which we present here, can be generalized to the field 
theory case. Therefore this toy model, which can be com- 
pared with a direct solution of the Schrodinger equation, 
is an ideal problem to test the accuracy of the numerical 
methods needed for field theory calculations. 



II. THE GENERATING FUNCTIONAL 

We consider a system of N oscillators of mass m cou- 
pled quadratically to a massless oscillator with coupling 
e. This quantum mechanical system is a model of a single 
momentum mode of scalar quantum electrodynamics [|S| . 
The Lagrangian for this system is given by: 
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We introduce the scaled variables, defined as 
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From now on, we use scaled variables. We wish to con- 
sider the time evolution of expectation values of observ- 
ables for an initial value problem, < t < oo. The way 
to formulate an initial value problem in quantum me- 
chanics, using a generating functional, was done more 
than thirty years ago by Schwinger, Bakshi and Mahan- 
thappa, and later by Keldysh j^]. This formalism, which 
is in the Heisenberg picture, is related to the fact that 
in the Schrodinger picture the evolution of the density 
matrix in quantum mechanics, 
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requires both a forward evolution from zero to t and a 
backward one from t to zero. The average value of ob- 
servables are given by traces over states of the system: 

{6(t)} = Tr{p(t)0} = Tr{p(0)O(t)} . 



in the Schrodinger and Heisenberg picture, respectively. 

As explained in the Appendix of ref. M, this neces- 
sitates both positive and negative time ordered opera- 
tors in the evolution of the observable operators and the 
introduction of two currents into the path integral for 
the generating functional. A concise way of writing the 
needed functions is to define the time integrals along a 
path in the complex time plane. This closed time path 
(CTP) is shown in Fig. [l]. The CTP integration contour 
is given by: 



F(t) dt 
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F-(t)dt. (2.4) 
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Using the CTP contour, the generating functional for the 
causal Gree n's f unctions for the theory described by the 
Lagrangian (2.1) is given by the path integral: 



Z[J,j] 



d[A] d[<p a ]e iNS ^ J ^ 



S[A,<t>;J,j}= JdtL, 
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The full closed time path Green's function for the two 
point functions is: 

G ab (t,t') = G> b (t,t')e c (t,t') +G<(t,t')Qc(t',t) , 

in terms of the Wightman functions, 

G> b (t,t>) = i{fa.(t)&(f)> - (Mt))(Mt'))} 
G< b (t,t') = t{(Mt')Mt)) - (Mt'))(Mt))} , 

where (0 o (t)&(f)) = Tr{p(0) 0«(t)<fc,(f )}> and where 
the CTP step function 0e(i, f) is defined by: 
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This is equivalent to a 2 x 2 matrix Green's function on 
the vector space {+, — }, often found in the literature. 

The large N expansion is obtained by performing the 
Gaussian integral over the cj> a variables to obtain an ef- 
fective action, and evaluating the remaining integral over 
A by the method of steepest descent. This gives: 
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where Ses[A ; J J] is given by 
S e ff[Ao; J, j] = 
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The stationary point Ao(t) is determined by 

( ,9 N 



di 2 
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A (t) = J(t). 

(2.9) 



Here <f>$ is defined as: 



— + (m 2 + e 2 Al) \ 0o o (t) - j„(t) , (2.10) 



and Go is given by the solution to 
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4>0a(t) and Ao(£) are to be regarded as functionals of 
the sources, J(t) and j(t). We have defined Sc(t,t') — 
dQ c (t,t')/dt. 

The inverse propagator Z) _1 (t, i') is 
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where 
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We solve ( |2.11 ) by introducing a complete set of solu- 
tions f(t) to the homogeneous equation, satisfying the 
Wronskian condition 

/«*(*)/-(*)-/:(*)/»(*) = -<■ (2.i5) 

The causal Green's functions Go can then be written: 

Goab(t,t') = 

iSab{fa(t)fZ(t')ec(t,t') + f a (t')fZ(t)6c(t',t)} . (2.16) 
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It is the solution of this equation that we will compare 
with our large N equations for the expectation values. 
Since at leading order in large N an initial Gaussian wave 
packet stays Gaussian under time evolution, we will start 
our problem with Gaussian initial data. Also we need to 
relate the parameters of the wave function at time zero 
to the values of one and two point functions and their 
time derivatives, since the Green's functions obey second 
order differential equations. At t — we choose our wave 
function to be a product of Gaussians: 
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where 
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The variables have the following meaning: 
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PA (o) = (-id/dA) t =o = (i(o)) 

we will evolve these equations (for N = 1) using a sym- 
plectic integrator as described in reference to compare 
with the results of the large N expansion. The param- 
eter G(t) in the Schrodinger wave function which is the 
real part of the width of the wave function is related to 
Go(t,t') of the Heisenberg approach via G(t) = Go(t,i)/i. 



III. TIME EVOLUTION EQUATIONS 

In the Schrodinger picture the Schrodinger equation 
rovcrns the time evolution of the N + 1 oscillators : 



A. Leading order in 1/N 

Let us now consider the equations and the initial con- 
ditions for the large N expansion. For simplicity and also 
because we are considering these equations as the single 
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mode approximation to scalar electrodynamics we con- 
sider the case where 4>o a = for all t. In terms of the 
mode functions f a discussed above we then have: 

G 0ab (t,t)/i =S ab f a (t)f:(t) (3.3) 

The coupled equations which need to be solved are: 
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The four initial conditions that have to be specified in 
lowest order are 

A(0); i(0); G(0), G(0) . 

Here we have assumed that the A field can be treated 
classically, which is the limit where the width of the fluc- 
tuations in A, namely D 7 can be ignored. In general we 
also have to specify D(0) and D(0). In terms of the f a 
this translates into the initial conditions 
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for all a. 



B. Order 1/N corrections 

The 1/N corrections to the generating fun ctio nal Z 
require the evaluation of the second term in (2/7). We 
invert ( 2.1 2| ) and find: 



D(t,t') = D (t,t') 



dh dt 2 D (tM)^o{tiM)D{t 2l t') . (3.6) 
c Jc 



Since we have chosen 4>o a (t) — 0, we can solve Do(t,t') 
as follows. Find the two linearly independent solutions g 
and g* to the homogeneous equation: 
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satisfying the Wronskian condition: 

9*(t)g(t)-g*(t)g(t) = -i , (3.8) 
In terms of these solutions we have: 
D (t,t') = i{g(t)g*(t')ec(t,t') + g(t')g*(t)G c (t',t)} 

(3.9) 



The initial width of the wave function is then given by 

D(0) = £> (0,0)/i= | 5 (0)| 2 (3.10) 

Thus we can relate the initial conditions on g to the initial 
conditions D(0) and D(0) of the wave function as follows: 
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a=l 

(3.11) 

Our numerical strategy for solving eq. (3.6) is discussed 
in appendix A. 

There are two ways to calculate the order 1 /N correc- 
tions to the time evolution problem. The first way is by 
a straight forward perturbation expansion of W[J, j] in 
powers of 1/N. If we consider the average value of A(t), 
we have: 
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Computing the derivatives, we obtain: 
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where 
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Here, all functions are to be evaluated using the first 
order solutions, Ao(t). 



C. Effective Action approach 

A second method of evaluation of the order 1/N cor- 
rections is to use a Legendre transformation to find the 
effective action F[A, 4>k\ 

T[A, 0] - W[J,j] - jf dt | J{t)A{t) + J2j a (t)<p a (t)^ , 

(3-14) 

Using the fact that S'cffL^; J, j] is a stationary point, we 
find, to order 1/A, 
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which is the classical expression plus the trace-log terms. 
Here, G~£[A](t,t') and D~ l [A, 0](i, t') are functional of 
the full A(t) and 4> a {t). Again, in this case we set 4> a {t) = 
0. This effective action agrees up to order 1/N with 
the true Legendre transform of the generating functional. 
One can now ignore where this action came from and 
directly write non-perturbative equations for A directly 
from this action which will agree to order 1/N with the 
previous equations. However at N = 1 where we will be 
making a numerical comparison with the exact answer, 
the results are expected to be quite different. We would 
like to know how the two results for A differ and which 
is more accurate at modest N. 

The equation of motion for A(t) which follows from 
varying the effective action T is: 
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The equation for A has to be solved simultaneously with 
the equations for G and D which now depend on the full 
A. To obtain G we now need to solve for the modes / 
which satisfy: 
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G is again given by: 



G a b(t,t') — 

lSab{fa(t)f:(t')e C (t,t') + f a (t')f* a (t)ec(t',t)} . (3.21) 

In terms of the full A and t hese new / we determine D 
using the integral equation (3.6) where now Dq is deter- 
mined using the new / in the equation for g. 

These equations of motion agree with the previous ones 
to order 1/N 2 . We had hoped that these equations which 
partially resum 1/N corrections would be more accurate 
at late times. However this turned out not to be the case, 
and in fact since A becomes unbounded sooner in this 
second approach, the late time behavior of this approach 
is worse. At N = 1 both methods agree with the exact 
answer for the same amount of time. When they diverge 
from each other they also diverge from the exact answer. 



IV. ENERGY 



The expectation value of the energy is given by: 
1 i N 

E/N = I(i 2 ) + i^{(0 a 2 ) +TO 2 (^ 2 ) + e 2 (A 2 ^)} 
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The expectation values are related to the full connected 
Green's Functions G,D K^ b , G A ah by the following equa- 
tions: 
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We notice that to order 1/N the connected 3 point 
function K 3 contributes but not the connected 4 point 
function G 4 . If we are using the generating functional 
approach we can now directly expand the energy in a 
power series in 1/N by assuming G n = Gq + -^G" + ... 
We then obtain: 

(A 2 (t)) = Al(t) + ^AoWMt) + jN D{t ^ 
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If instead we are using the effective action T to deter- 
mine the equations of motion, then we find the following 
expressions for the expectation values of the operators: 
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Where at order 1/N 
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(4.10) 
(4.11) 
(4.12) 
(4.13) 

G aa (t,t)D(t,t) 
(4.14) 



where now we use the full A and not Aq in determining 
all the Green's functions in K 3 . The energy for the T- 
method can then be calculated as 
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The above lead to the following expression of the energy: 



E z (t) = E (t) + -i^(t) 
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For Eo, we find 
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Using the equations of motion, we c an show t hat all ex- 
pressions for the energy, Eqs. 4/7, 4.15| ), are time 
independent, so that the energy is conserved for both 
the perturbative calculation as well as the effective ac- 
tion method. 



V. NUMERICAL RESULTS 

In order to find numerical solutions to these equations, 
we expand all functions in Chebyshev polynomials, and 
solve the resulting finite set of equations. In most of our 
calculations, we used a set of 32 polynomials. Details of 
this calculation are given in Appendix B. 

In order to make contact with our previous lowest or- 
der N = oo results we will use the set of initial conditions 



G 



(e = 0.3, E = 1, m = 1) given in reference 0. Time is 
measured in units of 1/m. We would first like to compare 
the two methods at large but finite N starting off with 
Gaussian initial data for the wave function. In a mean 
field approximation, such as lowest order large N, a wave 
function which starts out as a Gaussian remains Gaus- 
sian. Thus in mean field theory the only physical measur- 
able quantities are (A(t)), ((A(t) - (A(t))f) = D(t,t)/i, 
and (4> 2 (t)) = Q(t,t)/i. It is just these three moments 
that we will plot here to make contact with mean field 
results. (At infinite N, D is also zero, but it is finite in a 
Hartree approximation, as discussed in 0). We will de- 
note the straightforward use of the l/N expansion from 
the path integral as the Z-method and the results coming 
from the effective action approach, the T-method. 

In Fig. § we plot (A(t)) for N = 100, 16, and 8. We no- 
tice that as we reduce N the two different methods start 
diverging after one period. We also find that the second 
method which is based on the effective action has correc- 
tions which become unbounded at late times which leads 
us to believe that the second method is less stable. Both 
methods seem to have secular behavior (i.e. corrections 
which grow as a power in time). 

In Fig. ||, we show results for (4> 2 (t)) — Q(t,t)/i for 
N = 100, 16, and 8. We notice that at late times, this can 
go negative. This is a result of the fact that the operator 
<j> has a l/N expansion <f» = <f>a + jr(f>i + • • •• So that the 
positive definite width of the wave packet is ((^o + ^^i) 2 ) 
which also includes l/N 2 corrections. So once the l/N 
correction becomes as large as the lowest order term, the 
expansion breaks down and keeping terms only to order 
l/N in the expansion of the expectation value can lead 
to negative results. This negative result occurs after the 
two methods diverge which is an indication of when the 
expansion is breaking down. 

Now let us look at the effective width of the wave func- 
tion of the A oscillator. For the Z-method, we ha ve t hat 
the correlation function D(t,t'), defined by Eq. ( |3.6| ), is 
independent of N. However, the width determined in 
the effective action approach has l/N corrections due to 
the implicit dependence of A(t) on N. First, in Fig. 0, 
we compare the two approaches for calculating D(t, t)/i 
at N = 100. In Fig. ||, we show the N dependence of 
D(t,t)/i using the T-method at N = 100,16, and 8. 
From Fig. || we see that when t > 5 the l/N 2 effects 
are starting to appear. 

For the case N = 1, it is possible to compare the re- 
sults of both expansions with an exact numerical sim- 
ulation of the Schrodinger equation which we obtained 
in reference JtJ. Fig. ^ shows the time evolution of A(t) 
as computed using both the T- and the Z-method, com- 
pared with the exact solution for < t < 25. For com- 
parison we also include the lowest order in large N result. 
Both methods agree with each other and with the exact 
result for t < 8. They are accurate for at least twice the 
time scale as Ao(t), the first order large N result. It is 
gratifying to see that the two solutions diverge from each 



other approximately the same time as when they diverge 
from the exact solution so that the divergence of the two 
approaches which signals the onset of l/N 2 corrections 
sets the time scale for the accuracy of the result even 
for N = 1. 

In Fig. fj], we plot the effective width of the <j) oscillators 
((4>(t) — (<fi(t))) 2 ) = Gaa(t,t)/i using both methods and 
compare these results to the exact solution for N = 1 as 
well as the lowest order result. Here we notice a marked 
improvement of the l/N corrected results from the lowest 
order ones. We find again that the time at which the 
two methods start diverging from one another they also 
diverge from the exact answer. 

In Fig. |[ we plot the effective width of the A oscillator 
D(t,t)/i for the Z- and T-methods, as a function of time. 
In this case the lowest order in l/N is equivalent to a 
delta function width and so did not give a prediction. 
Here again we find the two solutions diverge from one 
another just when the approximation breaks down. 

What we can conclude from the above results are that 
l/N 2 corrections become important rather early in this 
particular time evolution problem. The period of agree- 
ment of our two methods is approximately the period 
when they give accurate results for the time evolution. 
We find at N = 1 that the width of the cj) oscillator is 
much better described when we add the l/N corrections. 
Also the time evolution of the A oscillator is now accu- 
rate for about twice the time period found previously at 
lowest order. 

In our calculations, we were able to verify energy con- 
servation to one part in 10 4 . In spite of this, when 
the l/N expansion breaks down it fails to preserve the 
posit ivit y of certain expectation values. If we look at 
Eq. J3.6| ) for D(t,t)/i, which is the positive definite ex- 
pectation value ({A{t) — (A(t)}) 2 ), we notice that this 
is an integral equation which sums all the bubbles and 
thus correctly takes into account the shift in the fre- 
quency of the quantum fluctuations of the A oscil lator. 
However in the equation for (A(t)) itself, Eq. ( |3.I6| ), the 
term which is the time dependent mass of the oscillator, 
namely (<p 2 (t)) , whi ch should also be positi ve defi nite, is 
given by Eq. ( fTj] ). We notice that Eq. ( |3.I7| ) is not 
an integral equation but only the leading term in a 1 /N 
expanded Green's function. This quantity need not be 
(and is not) positive definite. We can understand how 
this can happen in an example. In the vacuum state 
(<j> 2 (t)) = G(t,t)/i is time independent and is positive 
definite. One of the effects of the interactions is to renor- 
malize the mass (change the frequency) of the oscillator. 
So consider the toy problem: 



G(t,t) 



1 

2^ 

1 

2 



(uj 2 + ', 



N 



1-0 T - 

-1/2 



where we have kept only the e 2 /N correction to the mass 
shift of the cf> oscillator. Now m\ is positive since we have 
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assumed the interaction is repulsive. Again reexpanding 
G we get 



G(t,t) 



1 



2mo 



1 



2Nml 



So we see that if is large enough, then this quan- 
tity which is positive definite can appear negative when 
reexpanded in 1/N. 

This problem can also be a source of secular terms. 
For example consider the fact that the frequency of the 
classical A oscillator is shifted by 1/N corrections so that 
one has 



A(t) = A(0) cos (w t + ^rtj 



(5.1) 



If however the actual 1/N expansion for the oscillator 
instead gave the first two terms in the Taylor series we 
would get: 



A(t) = A(Q)[cos(w t) 



~N~ 



sin(a;o£)] 



(5.2) 



Using this expression would lead to blow up of A{t) in a 
time scale N/w\. This type of behavior is found in our 
numerical simulations. One way of solving the secular 
problem as well as guaranteeing the positivity of expecta- 
tion values is by ap propriately summing the series. That 
is we replace ( 3.17| ) by: 



Gab(t,t') = G a b(t,t) 

J2 [ dt i l^hG ac {tM)Vcd{tiM)Gdb(t 2 ,t') . 

(5.3) 



i 

N 



c Jc 



We th en ha ve to modify appropriately Eqs. ( 2.13 ), ( [2.14 ) 



( 3.18 ) and ( 3.19 ) in order to guarantee energy conserva- 
tion again. 

These modifications only change the results at order 
1/N 2 but are obviously very important because of the 
early breakdown of the naive methods. However in this 
paper we did not want to go beyond what one obtains di- 
rectly using the two direct approaches to the 1/N expan- 
sion and will discuss the resummed theory in a separate 
paper. 

This problem with expectation values becoming un- 
bounded is more acute in the T method. In Fig. [| we 
show that energy is conserved. However in Figs. [Ujand 
[TT] we show the above phenomena that several of the in- 
dividual contributions to the energy become negative if 
we keep terms only up to 1/N. 

The appearance of secular terms in a perturbation ex- 
pansion is quite well known jlC| ]. As an example for the 
classical anharmonic oscillator 



d 2 y 
dt 2 



-nf + v + gy = o , 



v = yo(t) +gyi(t) , 

and perform a perturbation series in the coupling con- 
stant g, secular terms linear in t appear in yi(t). One 
can see by doing a large- N expansion in the Schrodinger 
picture, that one can avoid the secular problem in a sim- 
ple way. If we solve the Schrodinger equation at large N 
in terms of the eigenfunctions and eigenvalues we find: 



A,t)=^ Cm^m((f>, A)e 



iE m t 



(5.4) 



if we assume a solution of the form 



where the ^ m ((j),A) are the solutions of the time inde- 
pendent Schrodinger equation. We can do a large N ex- 
pansion for the eigenfunctions and eigenvalues separately. 
We then find that ^f m (</>, A) has a Taylor series in 1/ \HV 
and E m has a Taylor series in 1/N. As long as we do not 
reexpand the exponentials having the time dependence 
(i.e. keep exp{— ilE^t + jjE^t ■ ■ ■]} ) then there are no 
secular terms. However if we reexpand the exponentials 
and keep only terms of order 1/N, secular terms propor- 
tional to powers of t will appear. We believe that this 
might very well be occuring in the large N expansion 
based on the CTP formalism. 



VI. CONCLUSIONS 

We have presented two methods for calculating the 
1/N corrections to the time evolution of a system of N+l 
oscillators using Schwinger's CTP formalism. These two 
methods differ by terms of order 1 /N 2 so that the diver- 
gence of results is an indication of the time during which 
keeping terms up to order 1 /N is accurate. This was veri- 
fied by comparing with direct numerical simulation of the 
Schrodinger equation for the case N = 1. We found that 
for N = 1, keeping the 1/N corrections allows us to ex- 
tend the range of time for which the expectations values 
track the exact answer significantly. This was most not- 
icablc for the the quantity which describes the width of 
the 4> wave function. 

In performing these numerical simulations, we found 
some shortcomings in the 1/N expansion. It does not 
guarantee the boundedness of various expectation val- 
ues, even though energy is exactly conserved. As a result 
of this, we expect that we are seeing secular behavior in 
our expansion which needs to be cured. One way of cur- 
ing this problem is to instead look at a 1/N expansion 
for the eigenfunctions and eigenvalues of the Schrodinger 
equation instead of working with the CTP formalism. 
Another is to resum the Green's function for the <f> os- 
cillator and modify the equations for the other variables 
appropriately so that energy is conserved. These two 
ideas will be discussed in a future work. 

In conclusion, by comparing our two direct approaches 
to performing a 1/N expansion, one using the generat- 
ing functional Z and one the effective action T, we found 
that the direct perturbation theory in 1/N from Z gave 
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results that were bounded for longer times. However hav- 
ing both methods, allowed us to determine when 1/N 2 
corrections became important. Thus having both meth- 
ods will be quite useful in simulating field theories, when 
there will be no exact solution to compare with. Solv- 
ing the secular problem wc found here will be discussed 
elsewhere. 
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APPENDIX A: EQUATION FOR D(T,T') 

In this appendix we discuss how we approach the prob- 
lem of finding the solution of eq. (|3.6|) 



Q>,<(t,t') 



Q*{t,t') + Q <t> (t',t). (A8) 



D(t,t') = D (t,t') 

- [ dh [ dt 2 D (t,t l )U(t 1 ,t 2 )D(t 2 ,t'). (Al) 
Jc Jc 

We first note that the causal Green's functions are sym- 
metric in the sense that A > {t,t') = _4<(t',i), and obey 
the additional condition 

■/*>,<(*, = -A* <t> (t,t') = A <t> (t',t). (A2) 

The last equation gives 

ne{A> (t, t')} = - Ke{A< (t, t')} (A3) 



lm{A > (t,t')} = - Jm{^<(t,t')} > 



or 



A>(t,t') - A^fat') 
A>(t,t') + A%(t,1?) 



2Tle{A > (t,t')} 
21m{A > {t,t')} . 



(A4) 

(A5) 
(A6) 



With these equations in mind, the causal Green's fum- 
ction A(t,t') is fully determined by the component 
A>(t,t') = 'R,e{A > (t,t')}+ilm{A > {t,t')}, and we need 
to evaluate the function only for t 1 < t. 

Directing our attention now to the calculation of 
D(t,t'), it is convenient to introduce 



Q(t,t') = / dt"D (t,t")W,t') 



(A7) 



Note that even though the functions Do(t, t 1 ) and II(i, t 1 ) 
satisfy the properties of the Green's functions listed 
above, the new function Q(t,t') satisfies only some, i.e. 



This does not prevent the Green's function D(t,t') to 
recover all the desired properties after one more CTP 
integration. In practice, this translates into the fact that 
we have to actually calculate the function Q > (t,t') for 
all moments t and t' , whereas the other causal functions 
need only so me h alf of the data. 
Then, eq. (Al) becomes 



D(t,t') = D (t,0 - J dt"Q{t,t")D(t",t'), (A9) 



D>(f,f) = D 0> {t,t') 

dt" {Q>(t,t")-Q<(tX)} !>>(*", *0 



+ / dt"Q>(t,i") [!>>(*", t')-D < (t",t')\ (A10) 
Jo 

£»<(*,*') = D <(*,*') 



dt" [Q>(t,t")-Q< (*.<")] D^f'X) 



+ / di"Q < (M")p>(i",t')-5 < (i",i')]. (All) 
Jo 



Using the properties (A2) - (AC), w e see that eq. ( All ) 
is redundant, and we can write eq. (A10) as 



£>>(M') = D 0> (t,t') 

-2 [ dt"Ke{Q > (t,t")} D>(t",lf) 
Jo 

+ 2 f dt'QyltX) ^{D>(t" \t')} . (A12) 
Jo 

We separate now the real and the imaginary part of (|Al| ) 
and obtain the system of integral equations 

UeiD^t')} = TZe{D 0> (t,t')} 

-2 [ dt'"R,e{Q > (t,t")}ne{D > (t",t')} 
Jo 

+ 2 f dt"Ke{Q > (t,t")}Ke{D > (t",t')} (A13) 

ImiDyit^')} = lm{D Q> (t,t')} 

-2 f dt"Ke{Q > (t,t")}lm{D > (t",t')} 
Jo 

+ 2 f dt"lm{Q > (t,t")}TZe{D > (t",t')} . (A14) 
Jo 

The above system of equations has to be solved for t' < t. 
Notice that the two equatio ns are independent, which 
allows us to solve first eq. (A13) for the real part of 
D > {t 1 and then use this resul t to d erive the imaginary 
part of D > (t,t') from equation ( A14). 
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APPENDIX B: NUMERICAL METHODS 

Our numerical technique involves the expansion of all 
the unknown functions, f(t), g(t), A(t), D(t,t'), in a 
Chebyshev polynomial basis. We follow a method devel- 
oped by El-gendy JO) and have applied it to combined 
differential and integrals equations of the type we have 
here. We use the same Chebyshev expansion methods 
for solving the Green's function equation for D(t,t'), in 
the CTP formalism, as explained in appendix A. In ad- 
dition, we divide the time up into small blocks, moving 
along block by block. 

Chebyshev polynomials of the first kind of degree n are 
defined by, 



T n {x) = cos(n arccos x) 



(Bl) 



We define the grid in the interval [—1,1] by choosing the 
positions: 



7 7T 

cos — 7 
N J 



0,1,. 



(B2) 



Then, the Chebyshev polynomials satisfy a discrete or- 
thogonality relation of the form 



N 

E 

fc=0 



Ti(x k )Tj(x k ) = Pi Si 



(B3) 



where the constants Pi are 
N 

Pi 



2 

N, 



i = 0,N 



(B4) 



An arbitrary function f(x) can be approximated in the 
interval [-1,1] by the formula 



N 

E 

3=0 



" a j T j( x ) , 



(B5) 



We denote f k = f(xk)- The coefficients dj are defined by 



2 N 

- X^'/OWfa) , 3=0,. 

k=Q 



,N (B6) 



and the summation symbol with double primes denotes 
a sum with first and last terms halved. As shown by 
El-gendy, we never have to actually compute these ex- 
pansion coefficients aj. Instead, we find fj directly. The 
adva ntage of the Chebyshev expansion method is that 
(B5), which is an expansion in a finite set, is exact for 
x = Xj,j = 0,...,N. 

We can approximate the calculation of the indefinite 
integral f_ x f(t) dt by 



I(x) = 



f f(t)dt = £ "a, f Tj(t)dt. (B7) 

J - 1 j=0 ^ 



where the integral Tj (i) dt is given by 

( Tj + ±(x) T^jx) 

2(j + 1) 2(j - 1) + f - 1 3 ~ 

\ Mx) - 1] if j = 1 

lTi(a:) + l if j = 



(B8) 



When the variable x takes the values of the grid points 
(B2), we can write eq. (B7) in a matrix format, like 



N 



L = I(xi 



E" », x L- 

3=0 



(B9) 



where S' is a square matrix of order (N+l), whose 
elements -E^" 1 ' 1 ' are equal to 

^Tjv^T/v-iOkj) + 2N ^ + T N+1 (xi)T N (xj) 
N 2 (-) k+1 1 



fcO(Ml) 

1 



2N 



+ N(N- 1 ) r -^-i(^)[ T iv-2(^) - ^T N (x 3 )] 



jV-2 



E J^ T k(^)\ T k-l{xj) - T k+1 (Xj)] 



k=l 



(BIO) 



In a similar manner, the derivative df(x)/dx is approx- 
imated by 



D(x) 



d 

dx 



N 



f{ X ) = J2" a 3 
3=0 



(Bll) 



where T'j(x) is the derivative of the Chebyshev polyno- 
mial of the first kind of degree j. By inserting the expres- 
sion of the a,j coefficients, and evaluating the derivative 
at x — it we obtain a matrix equation, 



N 



D, 



D{ Xl ) = 



// 5[-l.l] 



13 ■ 



(B12) 



3=0 



Here is given by: 



N 



51-1,11 = |fE" t '^) t ^) 



D 



(B13) 



fe=0 



Note that the approximations ( JBSj ) and (B12) are exact 
for x = = 0,...,N. 

Finally, the approximations (Bq, Bj| B12) can be gen- 
eralized by allowing the range of the approximation to 
be between two arbitrary limits a and b, instead of just 
-1 to 1. This is done performing the change of variable 
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(b-a) 



As a consequence, the matrices B and B become 
B [a,b] = b ~ a B [-i,i] 



g[a,b] 



b — a 



(B14) 

(B15) 
(B16) 



and the Chebyshev polynomials in eq. (B5) become 
now functions of the variable y. The matrices [/], 

f , f(t)dt and [df(x)/dx] will give then the value of 

the function f(x), its integral and derivative, at the co- 
ordinates Xk = Dkip — o)/2 + [b + a)/2, where yk are the 



(N+l) coordinates given by eq. (B2) 



t - plane 



e. 



FIG. 1. Complex time contour C for the closed time path 
integrals. 
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(a) N = 100. 
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(b) N = 16. 
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(c) N = 8. 

FIG. 2. Time evolution of A(t) for N = 8, 16, and 100, 
using the Z-method and the r-method. 




(c) N = 8. 

FIG. 3. Time evolution of G(t,t)/i for A/ = 
using the Z-method and the r-method. 



16, and 100, 
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FIG. 4. Time evolution of D(t,t)/i for N = 100, as com- FIG - 7 - Gaa(t,t)/i for the first and second order large N 

puted using the Z-method and the F-method. approximations as a function of time. 




FIG. 5. Time evolution of D{t,t)/i for N = 8, 16, 100, as 
computed using the F-method. 




FIG. 8. Comparison of exact numerical simulation of 
D(t, i) ji with the value computed using the Z and V methods. 




FIG. 6. Time evolution of A(t) for N = 1, as computed 
using the Z-method and the F-method, compared with the 
exact solution. 



n 



10 



15 



21) 



25 



t 



FIG. 9. The energy Ez = Eo + E\ /N, computed using the 
Z- method, and Er, using the F-method. 
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FIG. 10. Various components of the energy, as computed 

in the Z-method, as a function of time. Here we have in- 

1 1 

troduced the notations: E a = —A 2 (t), = - — 

1 



d 2 g aa (t,t') 

dtdt> 

i=t' 

1 d 2 D(t,t') 

2iN at9t ' 



t=t/ - ^ = S-i e2 <?„„(*,*) D(M), 

e 2 A(t) K 3 (t,t;t). 




FIG. 11. Various components of the energy, as computed 
in the Z-method, as a function of time. The notations are 
similar to those of Fig. llOl with the addition that now we 
have g aa (t,t') = GoooCMj + G Xaa {t,t')/N , and G aa (t,t') 
becomes Go aa (t,t'). 
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